function write_dynamics(accel)

%This function writes the dynamics file for the double pendulum.

filename = 'dynamics';

comments{1} = ['DZ = ' upper(filename) '(T,Z,U,P)'];
comments{2} = ' ';
comments{3} = 'FUNCTION:  ';
comments{4} = '    This function computes the dynamics for a point-mass that';
comments{5} = '    is connected to the origin by a force actuator.';
comments{6} = ' ';
comments{7}  = 'INPUTS: ';
comments{8}  = '    t = time. Dummy input for ode45. Not used.';
comments{9} = '    z = [4xn] matrix of states.';
comments{10} = '    P = struct of parameters';
comments{11} = 'OUTPUTS: ';
comments{12} = '    dz = [4xn] matrix of state derivatives';
comments{13} = ' ';
comments{14} = 'NOTES:';
comments{15} = ['    This file was automatically generated by ' mfilename]; 

params{1} = {'m','mass'};
params{2} = {'g ','gravity'};

states{1} = {'r','distance from point-mass to origin'};
states{2} = {'th','angle of vector between point-mass and origin'};
states{3} = {'dr','dr/dt'};
states{4} = {'dth','dth/dt'};

dstates{1} = {'dr','d/dt (distance from point-mass to origin)'};
dstates{2} = {'dth','d/dt (angle of vector between point-mass and origin)'};
dstates{3} = {'ddr','d/dt (dr/dt)'};
dstates{4} = {'ddth','d/dt (dth/dt)'};

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
%                               write file                                %
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~%
fid = fopen([filename '.m'],'w');

fprintf(fid, ['function dz = ' filename '(~,z,u,P) \n']);

for i=1:length(comments)
    fprintf(fid,['%%' comments{i} '\n']);
end
fprintf(fid,'\n');

for i=1:length(params)
    fprintf(fid,[params{i}{1} ' = P.' params{i}{1} '; %%' params{i}{2} '\n']);
end
fprintf(fid,'\n');

for i=1:length(states)
    fprintf(fid,[states{i}{1} ' = z(' num2str(i) ',:); %%' states{i}{2} '\n']);
end
fprintf(fid,'\n');

fprintf(fid,['ddr = ' vectorize(char(accel(1))) ';\n']);
fprintf(fid,['ddth = ' vectorize(char(accel(2))) ';\n']);
fprintf(fid,'\n');

fprintf(fid,'dz = [...\n');
for i=1:length(dstates)
    fprintf(fid,['    ' dstates{i}{1} '; %%' dstates{i}{2} '\n']);
end
fprintf(fid,'];\n');

fprintf(fid,'end \n');

fclose(fid);
end


